# This script makes all the plots based on GSS data (Figures 1, 2, 3, and C1 in the Appendix)
# Please note that we use dplyr to filter the dataset by plot, instead of making separate dfs

# Load packages and install them if necessary (automatic)
need = c(
  "dplyr",           #data structuring/cleaning
  "extrafont",       #plotting (fonts)
  "ggplot2",         #plotting
  "ggpubr",          #plotting (combining plots)
  "here",            #set wd
  "jtools",          #plotting (theme)
  "tidyr"            #data structuring/cleaning
)     
have = need %in% rownames(installed.packages())
if(any(!have)) {install.packages(need[!have])}
pack <- lapply(need, library, character.only = TRUE)
# Clean environment
rm(have, need, pack)

# Load data
load(file = "D:/dv_rep/Datasets/gss_boot_Plot.Rdata")

# Load fonts
extrafont::loadfonts()

#### FIGURE 1 (COHORT) ####

# Divergence
cohortDiv <- gss_Plot %>%
  filter(
    group == "Cohort" &
      boot == "Boot" &
      (facet %in% c("Divergence")) &
      (!(facet == "Sorting" & !grepl("^:", term)) | facet == "Sorting")) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("Greatest", "Silent", "Boomer", "Gen X", "Millennial"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  ylim(c(0.1, 0.27)) + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Deviation from predicted mean", x = "", col = "Controls: ", title = "Divergence") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), size = 3)),
    shape = "none",
    linetype = "none"
  )
# Sorting
cohortSort <- gss_Plot %>%
  filter(
    group == "Cohort" &
      boot == "Boot" &
      (facet %in% c("Sorting")) &
      facet == "Sorting" & grepl("^:", term)) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("Greatest", "Silent", "Boomer", "Gen X", "Millennial"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  scale_y_continuous(limits = c(-0.05, 0.60), breaks = seq(0, 0.60, by = 0.20), position = "right") + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Standardized effect of PID", x = "", col = "Controls: ", title = "Partisan Sorting") + 
  theme(legend.position = "bottom", 
        text = element_text(family = "Gill Sans MT"),
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 10, face = "bold"),
        legend.text = element_text(family = "Gill Sans MT"),
        legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )

# Extract common legend
legend_genMain <- get_legend(cohortDiv)
# Merge plots
gss_gen_Main <- ggarrange(
  cohortDiv + theme(legend.position = "none"),
  cohortSort + theme(legend.position = "none"),
  ncol = 2
)
# Merge with legend
gss_gen_Main <- ggarrange(
  gss_gen_Main,
  legend_genMain,
  ncol = 1,
  heights = c(20, 1)
)
# Annotate Figure
gss_Main_genPlot <- annotate_figure(gss_gen_Main, top = text_grob("Cohort effects", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))
# Save output (as .png and .pdf)
ggsave(gss_Main_genPlot, file = "fig1.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))
ggsave(gss_Main_genPlot, file = "fig1.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

#### FIGURE 2 (PERIOD) ####
# Divergence
yearDiv <- gss_Plot %>%
  filter(
    group == "Year" &
      boot == "Boot" &
      (facet %in% c("Divergence")) &
      (!(facet == "Sorting" & !grepl("^:", term)) | facet == "Sorting")) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term)) %>%
  ggplot(., aes(x = as.numeric(term), y = value, group = controls, col = controls, shape = controls, linetype = controls)) + 
  ylim(c(0.1, 0.27)) + 
  geom_point(position = position_dodge(width = 0.5), size = 2.5) +
  geom_line() + 
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  #scale_x_continuous(breaks = seq(1950, 2021, by = 10), labels = seq(1950, 2021, by = 10)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Deviation from predicted mean", x = "", col = "Controls: ", title = "Divergence") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), size = 3, linetype = c(1, 2))),
    shape = "none",
    linetype = "none"
  )
# Sorting
yearSort <- gss_Plot %>%
  filter(
    group == "Year" &
      boot == "Boot" &
      (facet %in% c("Sorting")) &
      facet == "Sorting" & grepl("^:", term)) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term)) %>%
  ggplot(., aes(x = as.numeric(term), y = value, group = controls, col = controls, shape = controls, linetype = controls)) + 
  scale_y_continuous(limits = c(-0.20, 0.80), breaks = seq(-0.20, 0.80, by = 0.20), position = "right") + 
  geom_point(position = position_dodge(width = 0.5), size = 2.5) +
  geom_line() + 
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Standardized effect of PID", x = "", col = "Controls: ", title = "Partisan Sorting") + 
  theme(legend.position = "bottom", 
        text = element_text(family = "Gill Sans MT"),
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 10, face = "bold"),
        legend.text = element_text(family = "Gill Sans MT"),
        legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), linetype = c(1, 2))),
    shape = "none",
    linetype = "none"
  )

# Extract common legend
legend_yearMain <- get_legend(yearDiv)
# Merge plots
gss_year_Main <- ggarrange(
  yearDiv + theme(legend.position = "none"),
  yearSort + theme(legend.position = "none"),
  ncol = 2
)
# Merge with legend
gss_year_Main <- ggarrange(
  gss_year_Main,
  legend_yearMain,
  ncol = 1,
  heights = c(20, 1)
)
# Annotate Figure
gss_Main_yearPlot <- annotate_figure(gss_year_Main, top = text_grob("Period effects", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))
# Save output (as .png and .pdf)
ggsave(gss_Main_yearPlot, file = "fig2.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))
ggsave(gss_Main_yearPlot, file = "fig2.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

#### FIGURE 3 (AGE) ####
# Divergence
ageDiv <- gss_Plot %>%
  filter(
    group == "Age" &
      boot == "Boot" &
      (facet %in% c("Divergence")) &
      (!(facet == "Sorting" & !grepl("^:", term)) | facet == "Sorting")) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("17-21", "22-29", "30-64", "65+"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  ylim(c(0.1, 0.27)) + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Deviation from predicted mean", x = "", col = "Controls: ", title = "Divergence") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17), size = 3)),
    shape = "none",
    linetype = "none"
  )
# Sorting
ageSort <- gss_Plot %>%
  filter(
    group == "Age" &
      boot == "Boot" &
      (facet %in% c("Sorting")) &
      facet == "Sorting" & grepl("^:", term)) %>%
  mutate(term = ifelse(facet == "Sorting", gsub("^:", "", term), term),
         term = factor(term, levels = c("17-21", "22-29", "30-64", "65+"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  scale_y_continuous(limits = c(-0.05, 0.60), breaks = seq(0, 0.60, by = 0.20), position = "right") + 
  geom_point(position = position_dodge(width = 0.5), size = 3.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0.2, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, scales = "free_y", repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "Standardized effect of PID", x = "", col = "Controls: ", title = "Partisan Sorting") + 
  theme(legend.position = "bottom", 
        text = element_text(family = "Gill Sans MT"),
        axis.text.y.right = element_text(hjust = 0.5), 
        plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(size = 8),
        axis.title.y = element_text(size = 10, face = "bold"),
        legend.text = element_text(family = "Gill Sans MT"),
        legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )
# Extract common legend
legend_ageMain <- get_legend(ageDiv)
# Merge plots
gss_age_Main <- ggarrange(
  ageDiv + theme(legend.position = "none"),
  ageSort + theme(legend.position = "none"),
  ncol = 2
)
# Merge with legend
gss_age_Main <- ggarrange(
  gss_age_Main,
  legend_ageMain,
  ncol = 1,
  heights = c(20, 1)
)
# Annotate Figure
gss_Main_agePlot <- annotate_figure(gss_age_Main, top = text_grob("Age effects", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))
# Save output (as .png and .pdf)
ggsave(gss_Main_agePlot, file = "fig3.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))
ggsave(gss_Main_agePlot, file = "fig3.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))

#### MEAN PLOTS (Appendix) ####
# Plot cohort effects
cohortMean <- gss_Plot %>%
  filter(group == "Cohort") %>%
  filter(boot == "Boot" & facet == "Mean" ) %>%
  mutate(term = factor(term, levels = c("Greatest", "Silent", "Boomer", "Gen X", "Millennial"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, ncol = 1, repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "<- Liberal | Conservative ->", x = "", col = "Controls: ", title = "Cohort") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )
# Plot year effects
yearMean <- gss_Plot %>%
  filter(group == "Year") %>%
  filter(boot == "Boot" & facet == "Mean" ) %>%
  ggplot(., aes(x = as.numeric(term), y = value, group = controls, col = controls, shape = controls)) + 
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_line() + 
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, ncol = 1, repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "<- Liberal | Conservative ->", x = "", col = "Controls: ", title = "Period") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )
# Plot age effects
ageMean <- gss_Plot %>%
  filter(group == "Age") %>%
  filter(boot == "Boot" & facet == "Mean" ) %>%
  mutate(term = factor(term, levels = c("17-21", "22-29", "30-64", "65+"))) %>%
  ggplot(., aes(x = term, y = value, group = controls, col = controls, shape = controls)) + 
  geom_point(position = position_dodge(width = 0.5), size = 1.5) +
  geom_errorbar(aes(ymin = lower_ci_value, ymax = upper_ci_value), width = 0, position = position_dodge(width = 0.5)) + 
  lemon::facet_rep_wrap(~ var, nrow = 5, ncol = 1, repeat.tick.labels = TRUE) + 
  theme_apa() +
  labs(y = "<- Liberal | Conservative ->", x = "", col = "Controls: ", title = "Age") + 
  theme(
    legend.position = "bottom", 
    text = element_text(family = "Gill Sans MT"),
    plot.title = element_text(hjust = 0.5),
    axis.text.x = element_text(size = 8),
    axis.title.y = element_text(size = 10, face = "bold"),
    legend.text = element_text(family = "Gill Sans MT"),
    legend.title = element_text(family = "Gill Sans MT")) +
  guides(
    col = guide_legend(override.aes = list(shape = c(16, 17))),
    shape = "none",
    linetype = "none"
  )
# Obtain common legend
legend_Mean <- get_legend(yearMean)
# Merge plots
mean_Plot <- ggarrange(
  cohortMean + theme(legend.position = "none"),
  yearMean   + theme(legend.position = "none"),
  ageMean    + theme(legend.position = "none"),
  ncol = 3
)
# Merge with legend
mean_Plot <- ggarrange(
  mean_Plot,
  legend_Mean,
  ncol = 1,
  heights = c(20, 1)
)
# Annotate Figure
mean_Plot <- annotate_figure(mean_Plot, top = text_grob("Mean-levels (GSS)", hjust = 0.5, face = "bold", family = "Gill Sans MT", size = 20))
# Save output as .png and .pdf
ggsave(mean_Plot, file = "figc1.png", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))
ggsave(mean_Plot, file = "figc1.pdf", dpi = 900, 
       width = 21,
       height = 29.7,
       units = c("cm"))a